R0Beta

input data

source(file = '../inputProcessing/date_input.r')


rep <- 1e2
# rep_sim <- 1e4


for (date_num in 1:length(dates_input)){
  
  date_week_finishing <-  as.Date(dates_input[date_num],format = '%Y-%m-%d')
  output <- list()
  
  # date_week_finishing <-  as.Date('29/03/2020',format = '%d/%m/%Y')
  # # si <- 2
  # rep <- 1e2
  # rep_sim <- 1e2
  
  for (Mdata in c('A','G')){
    
    inputs<- readRDS(file=paste0(
      '../../Mobility.2.0.Save/Saved_file/inputProcessing/Rdata/',
      date_week_finishing,'_inputs_from_',Mdata,'.rds'))
    
    for (si in 1:2){
      
      ########################################
      ####inputs
      D <- inputs$D
      M <- inputs$M
      Ot <- inputs$Ot[[si]]$Ot
      W <- inputs$W[[si]]$W
      H <- inputs$H
      SI <- inputs$SI[[si]]
      delta_id <- inputs$delta_id
      #
      country <- names(D)[-1]
      
      mD <- as.matrix(D[,-1])
      N_geo <- ncol(D)-1
      N_days <- nrow(D)
      
      N_week <- nrow(D)/7
      
      # prerp MCMC
      # epi
      R0 <- 5
      # for risk
      # b <- 0
      o <- 1
      
      theta0 <- c(rep(R0,N_geo*N_week),
                  rep(o,N_geo*N_week))#,
      # rep(b,N_geo))
      # theta0[2] <- 2
      # theta0[3] <- 3
      prior_theta <- matrix(c(rep(c(0,5),N_geo*N_week),rep(c(0,1e5),N_geo*N_week)),
                            length(theta0),2, byrow=TRUE)
      
      # parameter names
      f0 <- function(x) paste0('R',x)
      f1 <- function(x,y) paste0(x,'_',y)
      # f1 <- function(x) paste0('beta_',x)
      # f2 <- function(x) paste0('Over_',x)
      n_t<- rep(sapply(1:N_week,f0),N_geo)#, sapply(country,f1))#, sapply(country,f2))
      n_t<- rep(country,each = N_week)#, sapply(country,f1))#, sapply(country,f2))
      n_t<-c(n_t,'over')
      # sd dev for proposal
      sigma <- rep(1e-1,length(theta0))
      
      ########################################
      # useful functions
      sapply(paste0('Rscript/',(list.files('Rscript/'))),FUN = source)
      
      # R_mat <- matrix(rep(theta0,each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
      
      ########################################
      #run MCMC
      # res <- MCMC_iter(iter = rep, theta0 = theta0, s = sigma)
      res <- MCMC_full(iter = rep, 
                       theta0 = theta0,
                       s = sigma,
                       repli_adapt = 10,
                       within_iter = rep/10)
      
      
      
      ########################################
      ## check convergence
      
       print('########################################')
      print(paste0('check convergence',date_week_finishing,' - ',Mdata,' - ',si)) 

      
      Acc <- colSums(diff(res$theta)!=0)/(rep-1)
      # Acc
      if(rep >1e2) {
        f <- round(seq(1,rep,length.out = 1e2))
      }else{
        f <- 1:rep
      }
      plot(res$logL[,1],main=paste0('DIC ',res$DIC[1],', P ',res$DIC[2]))
      layout(matrix(1:4,2,2))
      for (i in 1:(length(theta0)/2+1)){
        if(i == (length(theta0)/2+1)){
          plot(res$theta[f,i],
               main = paste0(n_t[i],' - ',round(Acc[i]*100)) )
        }else{
          plot(res$theta[f,i],
               main = paste0(n_t[i],' - ',round(Acc[i]*100)) , ylim = prior_theta[i,])
          
        }
      }
      
      # apply(res$theta,2,quantile,c(.5,.025,.975))
      
      
      
      
      
      ########################################
      ## Baseline model: estimate Rt
      
      
      summary <- apply(res$theta,2,quantile,c(.5,.025,.975))
      
      median <- matrix(rep(summary[1,],each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
      low <- matrix(rep(summary[2,],each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
      up <- matrix(rep(summary[3,],each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
      
      mean_r <- matrix(rep(apply(res$theta,2,mean),each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
      sd_r <- matrix(rep(apply(res$theta,2,sd),each=7),nrow = N_days, ncol = N_geo,byrow = FALSE)
      
      # format outputs
      temp <- D
      temp[,-1] <- NA
      
      results_baseline <- list(median_R = temp,
                               low_R = temp,
                               up_R = temp,
                               CV = temp,
                               CV_below.2 = temp,
                               mean = temp)
      
      for (i in 1:N_geo){
        f <- which(cumsum(D[,i+1])>0)
        results_baseline$median_R[f,i+1] <- median[f,i]
        results_baseline$low_R[f,i+1] <- low[f,i]
        results_baseline$up_R[f,i+1] <- up[f,i]
        
        results_baseline$CV[f,i+1] <- sd_r[f,i]/mean_r[f,i]
        results_baseline$CV_below.2[f,i+1] <- results_baseline$CV[f,i+1] < .2
        results_baseline$mean[f,i+1] <- mean_r[f,i]
        
        # results_baseline$M_D[,i+1] <- M_D[,i]
      }
      results_baseline$DIC <- res$DIC
      
      
      output[[si]] <- list(resEpi = results_baseline)
      
    }
    
    
    
    saveRDS(object = output,
            file = paste0('../../Mobility.2.0.Save/Saved_file/EpiEstim/Rdata/',
                          date_week_finishing,'_output_from_',Mdata,'.rds' ))
    
  }
}
## [1] "########################################"
## [1] "check convergence2020-05-03 - A - 1"

## [1] "########################################"
## [1] "check convergence2020-05-03 - A - 2"

## [1] "########################################"
## [1] "check convergence2020-05-03 - G - 1"

## [1] "########################################"
## [1] "check convergence2020-05-03 - G - 2"

## [1] "########################################"
## [1] "check convergence2020-05-10 - A - 1"

## [1] "########################################"
## [1] "check convergence2020-05-10 - A - 2"

## [1] "########################################"
## [1] "check convergence2020-05-10 - G - 1"

## [1] "########################################"
## [1] "check convergence2020-05-10 - G - 2"